NASA  Contractor  Report  187542 
ICASE  Report  No.  91-32 


ICASE 


~r  ♦  n 


AD-A236  842 


THE  P  —  RKDG  METHOD  FOR  TWO-DIMENSIONAL 
EULER  EQUATIONS  OF  GAS  DYNAMICS 


Bernardo  Cockburn 
Chi-Wang  Shu 


Contract  No.  NAS  1-1 8605 
March  1991 


Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  Virginia  23665-5225 

Operated  by  the  Universities  Space  Research  Association 


I 

i 

NASA 

National  Aeronautics  and 
Space  Administration 

Langley  Research  Center 

Hampton,  Virginia  23665-5225 


/i 


i 


’.m 


91-01570 


91  6  7  09r 


THE  P 1  -  RKDG  METHOD  FOR  TWO-DIMENSIONAL 
EULER  EQUATIONS  OF  GAS  DYNAMICS 


Bernardo  Cockburn1 
School  of  Mathematics 
University  of  Minnesota 
Minneapolis,  MN  55455 
and 

Chi-Wang  Shu2 3 
Division  of  Applied  Mathematics 
Brown  University 
Providence,  Rhode  Island  02912 


ABSTRACT 


it. 

.  *  - 

^ - V 

■  L*T 

i  ^  <•;, 

\ 

■  Ji.’.! 

:  j.i 

;  sy 

Cl 

V-'*  uaul  •  •  V 

I  rnA/t’T 


| 

1 

cnee  led. 

lA 

i  I 

| 

We  continue  our  earlier  work  on  a  class  on  nonlinearly  stable  Runge-Kutta  local  projec¬ 
tion  discontinuous  Galerkin  (RKDG)  finite  element  methods  for  conservation  laws.  Two- 
dimensional  Euler  equations  for  gas  dynamics  are  solved  using  P1  elements.  We  discuss 
the  generalization  of  the  local  projection,  which  for  scalar  nonlinear  conservation  laws  was 
designed  to  satisfy  a  local  maximum  principle,  to  systems  of  conservation  laws  such  as  the 
Euler  equations  of  gas  dynamics  using  local  characteristic  decompositions.  Numerical  exam¬ 
ples  include  the  standard  regular  shock  reflection  problem,  the  forward  facing  step  problem 
and  the  double  Mach  reflection  problem.  These  preliminary  numerical  examples  are  chosen 
to  show  the  capacity  of  our  approach  to  obtain  nonlinearly  stable  results  comparable  with 
the  modern  nonoscillatory  finite  difference  methods.  Generalizations  to  Pfc  elements  with 
k  >  1  and  the  use  of  adaptive  triangulations  to  minimize  local  errors  constitute  ongoing 
research. 
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1.  Introduction.  In  this  paper  we  continue  our  earlier  work  [Cl],  [C2],  [C3],  [C4], 
of  constructing  and  analyzing  a  class  of  discontinuous  Galerkin  finite  element  method  for 
solving  conservation  laws 
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equipped  with  suitable  initial  or  initial-boundary  conditions.  We  concentrate  on  two- 
dimensional  Euler  equations  of  gas  dynamics,  i.e.,  in  (1.1)  with  d  =  2,  u  =  ( p,pqx,pqy,EY , 
fi(u)  =  gJ.u  +  (0,p,0,gIp)',  f2(u)  =  qyii  +  (Q,  Q,p,qy  p)‘,  where  qx,qy  are  the  velocity  com¬ 
ponents  in  the  x  and  y  directions,  p  is  the  density,  E  is  the  total  energy, 
p  —  ( 7  —  1)(£'  —  \p{q2x  4-  9y))  is  the  pressure,  and  7  =  1.4  for  the  air. 

One  distinctive  feature  of  our  approach  is  a  local  projection  limiting  which  borrows  the 
successful  non-oscillatory  finite  difference  methodology,  guarantees  total  variation  bound¬ 
edness  (TVB)  for  one-dimensional  nonlinear  scalar  equations  and  linear  systems,  and  yields 
a  local  maximum  principle  for  multi- dimensional  scalar  equations.  Another  feature  of  our 
approach  is  the  use  of  high-order  total  variation  diminishing  (TVD)  Runge-Kutta  type 
time  discretizations  [Si]  which  renders  the  scheme  explicit  (and  hence  fully  parallelizable) 
and  computationally  efficient.  The  general  framework  was. given  in  [C2]  for  the  nonlinear 
one-dimensional  scalar  case.  The  TVB  property  and  convergence  were  proven  for  general 
Pfc  elements  (which  give  rise  to  uniformly  (k  +  l)-th  order  accurate  schemes).  Numeri¬ 
cal  results  for  k  =  1  and  k  =  2  (second  and  third  order)  were  shown  which  gave  sharp, 
monotone  shock  transitions  and  uniform  high-order  accuracy  in  the  smooth  part  of  the 
solutions.  In  [C3],  we  applied  the  method  to  the  one-dimensional  Euler  equations  of  gas 
dynamics  by  designing  the  local  projection  in  the  local  characteristic  fields.  The  resulting 
scheme  was  proven  TVB  for  linear  systems  with  both  initial  and  initial-boundary  con¬ 
ditions.  Numerical  results  included  both  standard  shock  tube  problems  and  a  problem 
involving  the  interactions  between  a  Mach  3  shock  and  a  density  wave  -  a  prototype 
for  shock-turbulence  interactions.  In  all  cases,  we  obtained  results  comparable  to  those 
obtained  by  recent  nonoscillatory  finite  difference  methods  (e.g.,  [W],  [S2]).  The  crucial 
generalization  to  multi-space  dimensions  was  carried  out  in  [C4],  where  we  introduced  a 
local  projection  limiting  which  does  not  have  a  direct  counter-part  in  the  current  finite 
difference  schemes,  guarantees  a  local  maximum  principle  for  a  class  of  very  general  trian¬ 
gulation  (we  called  them  B-triangulations),  and  maintains  uniform  high-order  accuracy  in 
the  smooth  part  of  the  solution.  Numerical  results  showed  the  potential  of  the  scheme  in 
easy  handling  general  triangulations  and  boundary  conditions. 

The  organization  of  this  paper  is  as  follows.  In  §2  we  briefly  describe  the  formulation 
of  the  scheme.  Special  attention  is  paid  to  the  monotone  fluxes  (approximate  Riemann 
solvers)  across  the  edges  of  the  triangles  and  to  the  local  projections  associated  with  them, 
since  they  are  distinct  from  the  scalar  case  considered  in  [C4],  In  §3  we  present  numerical 
results  of  our  P*-RKDG  scheme  on  the  the  shock  reflection  problem,  the  forward  facing 
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step  problem,  and  the  double  Mach  reflection  problem.  We  end  with  some  concluding 
remarks  in  §4. 

2.  The  formulation  of  the  scheme.  Suppose  we  axe  solving  the  equation  (1.1) 
( d  =  2)  on  a  polygonal  domain  Q.  Let  Th  =  {K}  be  a  B-triangulation  of  f l  (see  [C4]  for 
the  definition),  and  set  V),  =  {p  G  L°°(fl)  :  p\x  G  P k(K)  VK  G  7/J  where  Pfc(iir)  is  the 
space  of  polynomials  of  total  degree  less  than  or  equal  to  k  on  K.  In  this  paper  we  only 
consider  k  =  1,  which  yields  second-order  accurate  schemes.  Notice  that  the  functions  in 
Vh  can  be  discontinuous  across  the  edges  of  the  triangles  in  Th-  For  scalar  equations,  the 
discontinuous  Galerkin  method  consists  in  finding  U/,  :  [0,  Tj  •— *■  Vj,  satisfying 


4  J  uh(t,x)vh(x)dx  +  I  f(uh(t,x))  •  vh(x)dT 

J  ./-at*  » 


(2.1) 


e&dK 


-  J  i(uh(t,  x))  ■  grad  Vh(x)  dx  =  0,  Vvh  €  Vh, 

K 


where  ne  K  is  the  outward  unit  normal  to  the  edge  e  and  f  =  (/i,/2).  This  is  obtained 
by  multiplying  (1.1)  with  a  test  function  t>/,  G  Vh,  integrating  over  a  triangle  K  G  Th, 
integrating  formally  by  parts,  and  replacing  u  by  its  approximation  it/,.  Since  Uh  can  be 
discontinuous  across  an  edge  e,  we  replace  f(uh(t,x ))  •  ne>x  by  a  monotone  flux 

(2.2) 


which  is  consistent:  hCifc(u,u)  =  f(u)  •  n e,K',  monotone:  he^(u,  v)  is  nondecreasing  in  u 
and  nonincreasing  in  u;  Lipschitz  continuous;  and  legitimate  as  a  flux: 

K, Uh(x‘*,( K>))  =  for  K'nK  =  e. 


The  line  integrals  in  (2.1)  are  replaced  by  the  two-point  Gauss  quadrature  rule,  and  the 
surface  integrals  replaced  by  the  mid-point  quadrature  rule.  The  resulting  O.D.E.  is  then 
discretized  by  the  second-order  accurate  TVD  Runge-Kutta  method  in  [Si]  and  coupled 
with  a  local  projection  applied  at  the  end  of  each  inner  step.  This  local  projection  is 
crucial  for  keeping  the  nonlinear  stability  (a  local  maximum  principle)  of  the  scheme;  it  is 
described  in  detail  in  [C4]. 

For  systems  of  equations,  (2.1)  is  satisfied  by  each  component  of  u*,.  The  monotone 
flux  (2.2),  however,  is  replaced  by  an  approximate  Riemann  solver:  the  relevant  quantities 
f (uh(xtrit(K\  t)  ■  ne  K  and  i(uh(xext{K) ,  t)  ■  ne  K  are  left-multiplied  by  the  left  eigenvectors 
of  some  average  Jacobian  across  e  (we  use  the  Roe  average  [Roe]  associated  with  u k  and 
\iK' ,  the  cell  averages  of  u  in  K  and  K'  whose  common  edge  is  e,  in  the  normal  direction 
to  e),  the  scalar  monotone  flux  (2.2)  (we  use  the  local  Lax-Friedrichs  flux  [C2])  is  then 
applied  in  each  of  the  resulting  local  characteristic  fields,  and  the  results  are  projected  back 


2 


by  using  the  right  eigenvectors  of  the  same  Jacobian.  This  is  essentially  a  one-dimensional 
characteristics  procedure  across  the  edge  e.  Details  about  it  can  be  found  in  [C3]. 

Next  we  describe  the  local  projection.  For  the  scalar  case,  a  local  maximum  principle 
cam  be  proven  if  we  compute  local  approximate  gradients  using  and  u/cs  where  K  amd 
K'  are  the  immediate  amd  some  further  away  neighbors  of  AT,  amd  use  them  to  limit  the 
deviations  of  u/,  evaluated  at  Gaussian  points  from  the  mean  u/<-.  The  limiting  is  performed 
in  such  a  way  as  not  to  destroy  the  formal  accuracy  in  smooth  regions.  The  details  amd 
proofs  can  be  found  in  [C4] .  For  systems  of  equations,  the  limiting  should  be  performed  in 
the  local  characteristic  fields.  We  adopt  the  simplified  version  of  the  local  projection  which 
involves  only  the  midpoints  of  the  edges  of  K  instead  of  aill  the  Gaussian  points  (as  in  [C4] ) 
amd  limits  the  deviation  by  considering  only  the  approximate  gradient  obtained  by  using 
u/c,  uk"'!  and  u^<2,  where  K'\  and  K'z  are  the  two  ‘forward’  facing  triangles  adjacent  to 
A';  notice  that  the  vector  M  —  B/^,  where  is  the  baxy center  of  the  triangle  A',  is  a 
positive  combination  of  the  vectors  B/cq  —  B/c  and  B/^j  —  B*-,  see  Fig.  1.  (In  [C4]  we 
used  a  projection  that  uses  an  additional  approximate  gradient  obtained  with  backward 
facing  elements  I\\  and  K '2.  Our  numerical  experience  shows  that  such  a  projection  is 
not  needed  in  this  framework.)  The  relaxation  factor  6,  which  is  the  ratio  allowed  of  the 
deviation  of  u/,  to  u/c  over  that  computed  by  the  approximate  gradient,  is  chosen  to  be  1.5 
(  6  >  1  guarantees  second-order  accuracy  in  smooth  monotone  regions).  This  simplified 
projection  reduces  the  computational  cost  and  seems  to  work  well  numerically.  We  again 
use  the  Roe  average  associated  with  0#  and  uk'h  where  KC\K'i  =  e,  along  the  normal  of 
e,  to  perform  the  characteristic  decomposition  for  the  projection  limiting  at  the  midpoint 
of  e;  see  Fig.  1.  The  projection  limiting  is  done  at  each  midpoint  independently  form  the 
projections  at  the  other  midpoints.  Then,  a  simple  readjustment  is  performed  to  enforce 
the  conservativity  of  the  method.  This  renders  the  projection  limiting  simple  and  efficient. 
In  all  our  computations  we  took  M  —  50;  see  [C4j. 

3.  The  numerical  examples.  We  use  three  standard  test  problems:  the  regular 
shock  reflection,  flow  past  a  forward  facing  step,  and  the  double  Mach  reflection  problem, 
to  display  the  behavior  of  our  P’-RKDG  scheme.  We  run  our  code  on  triangulations 
that  are  obtained,  essentially,  by  taking  the  usual  finite  difference  grid  and  dividing  its 
rectangles  into  two  triangles.  (In  principle,  this  should  produce  a  bigger  numerical  diffusion 
and  a  stronger  mesh  distorsion;  on  the  other  hand,  these  undesirable  effects  could  be 
counterbalanced  by  the  inherent  increase  of  degrees  of  freedom.  Indeed,  this  seems  to 
be  the  case.)  The  objective  is  to  verify  that  our  scheme  gives  results  comparable  to 
the  results  given  by  the  recent  nonoscillatory  finite  difference  schemes  for  these  standard 
test  problems.  Except  for  the  second  run  in  Example  1,  we  deliberately  do  not  use  the 
alignments  of  triangles  with  shock  structures  in  order  to  carry  out  a  fair  comparison  with 
finite  difference  schemes.  For  problems  involving  complicated  geometries  and/or  boundary 
conditions,  the  advantages  of  our  finite  element  schemes  over  finite  difference  methods 
would  be  more  significant.  Also,  for  problems  involving  both  shocks  and  complicated  flow 
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Fig.  1.  Projecting  uj, |/c fM.)  at  the  midpoint  M.  The  deviation  ujJ/e^Mj-UK 
is  compared  with,  b  times  the  deviation  at  M  obtained  by  the  approximate  gradient 
formed  vnth  iin,  u k\>  an ^  U/f'j  • 

structures,  P k  -RKD  G  schemes  would  show  their  advantage  in  high  resolutions.  Research 
along  these  directions  is  currently  being  carried  out. 

In  all  our  computations  we  have  used  the  CFL-condition  Afsup  A#  \e\/\K\  <  0.3, 
where  A/c  is  the  modulus  of  the  velocity  plus  the  sound  speed  evaluated  at  the  bary center 
of  the  traingle  K  We  use  the  CRAY2  at  the  Minnesota  Supercomputer  Center  to  carry 
out  our  computations.  The  code  is  carefully  written  so  that  most  of  the  major  loops  are 
fully  vectorized.  Due  to  the  local  structure  of  the  algorithm,  a  parallel  version  of  the 
CRAY2  code  could  be  written;  we  plan  to  do  this  in  the  future.  Our  graphics  has  been 
done  with  the  finite  element  package  MODULEF.  We  have  used  30  isovalues  in  all  our 
contour  figures. 

Example  1:  Regular  shock  reflection.  This  well-known  example  involves  an 
oblique  shock  reflecting  from  the  lower  boundary  of  a  rectangular  domain.  The  com¬ 
putational  domain  is  0  <  x  <  4.12829,0  <  y  <  1.  The  initial  condition  is  p  =  l.qx  = 
2.9,  qy  —  0,p  =  1/1.4  throughout  the  computational  domain.  The  boundary  conditions 
applied  are:  inflow  boundary  condition  at  x  =  0  (prescribe  all  components  of  u/,  with 
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values  same  as  the  initial  condition);  outflow  boundary  condition  at  x  =  4.12829  (all  com¬ 
ponents  axe  freely  flowing  out);  enforcing  the  post-shock  condition  at  the  upper  boundary 
y  =  1  (p  =  1.7,  qx  =  2.61932,  qy  =  — . 506339,  p  =  1.52824);  and  enforcing  a  Reflection 
boundary  condition  at  the  lower  boundary  y  =  0  (if  an  edge  e  of  a  triangle  is  at  the  lower 
boundary  y  =  0,  then  the  boundary  value  of  u/,  at  (x,  y  =  0)  is  chosen  to  be  the  same 
as  Uk  from  y  *-*  0+  for  p,  qx,  and  p  but  opposite  in  sign  for  qy.  The  exact  solution  to 
this  problem  is  an  incoming  shock  of  29  degrees  with  the  lower  boundary  and  a  reflected 
shock  of  23.28  degrees.  The  exact  solution  past  the  second  shock  should  be  p  =  2.68732, 
qx  =  2.40148,  qy  =  0  and  p  =  2.93413. 

In  Fig.  2b  we  show  the  pressure  contour  computed  with  a  refinement  of  the  triangu¬ 
lation  in  Fig.  2a  corresponding  to  the  uniform  Cartesian  grid  ‘Ax  =  Ay  =  1/40.’  Notice 
how  the  incident  shock  is  approximated  much  better  that  the  reflected  one.  This  is  due 
to  the  fact  that  the  triangles  are  partially  aligned  with  the  incident  shock  and  cut  the 
reflected  shock,  increasing  in  this  case  the  numerical  diffusion;  see  Fig.  2b.  In  order  to 
see  how  dramatically  the  result  can  improve  when  triangles  are  perfectly  aligned  with  the 
shocks,  we  show  in  Fig.  2d  the  pressure  contour  computed  with  the  triangulation  in  Fig. 
2c.  In  this  case,  the  L1 -error  (with  only  three  triangles)  is  only  ten  times  bigger  than  the 
L1 -error  produced  by  the  approximation  displayed  on  Fig.  2.b  (which  uses  800  triangles). 

This  test  problem  is  simple  in  structure:  three  constant  values  separated  by  two  shocks. 
We  use  it  to  test  (i)  the  nonoscillatory  property  of  our  simplified  local  projection,  (ii)  the 
behavior  of  the  numerical  reflecting  boundary  condition  at  the  lower  boundary,  and  (iii) 
the  effect  when  triangles  are  aligned  with  shocks. 


Fig.  2. a.  The  triangulation  ‘Ax  =  Ay  =  1/10’.  The  pressure  below  has  been 
computed  on  a  triangulation  four  times  finer  than  this  one. 
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Fig.  2.B.  Pressure  contour  computed  on  the  triangulation  ‘Ax  =  Ay  =  1/40’. 


Fig.  2.C.  The  triangles  of  this  triangulation  are  aligned  along  the  shocks  of 
the  exact  solution. 


Fig.  2.D.  Pressure  contours  computed  on  the  triangulation  of  Fig.  2.c. 


Example  2:  Flow  past  a  forward  facing  step.  This  is  one  of  the  two-dimensional 
problems  Woodward  and  Colella  [W]  used  to  test  the  behavior  of  various  finite  difference 

f  0.6,  if  y  <  .2 

schemes.  The  computational  domain  is  0  <  y  <  1,  0  <  x  <  <  .The 

l  3.0,  ify  >  0.2 

initial  condition  is  a  Mach  3  uniform  flow:  p  =  1.4,  qz  =3 ,qy  =  0,p  =  1  throughout  the 
computational  domain.  The  boundary  conditions  applied  are:  inflow  at  x  =  0;  outflow 
at  x  =  3;  and  reflecting  at  the  walls  y  =  1,  y  =  0,  x  =  0.6  and  y  =  0.2.  The  comer 
is  a  singularity.  In  [W],  Woodward  and  Colella  suggested  a  way  to  numerically  treat  the 
singularity.  We  display  results  without,  Fig.  3.b,  and  with,  Fig.  3.c,  such  a  treatment. 
Figs.  3.b  and  3.c  show  the  density  contours  at  T  =  4  computed  with  a  refinement  of  the 
triangulation  shown  in  Fig.  3. a  corresponding  to  the  Cartesian  grid  ‘Ax  =  Ay  =  1/40’. 
We  can  see  that  our  P1  -RKDG  scheme  produces  results  comparable  with  the  same  order 
MUSCL  finite  difference  scheme  using  the  grid  !Ax  =  Ay  =  1/40’..  We  remark  that  since  no 
sharpening  technique  is  applied  in  the  linearly  degenerate  field,  the  contact  discontinuities 
are  more  seriously  smeared  than  shocks.  We  are  currently  investigating  the  application  of 
Yang’s  artificial  compression  ideas  [Y]. 

Example  3:  Double  Mach  reflection.  This  is  the  second  example  used  in  [W]  to 
compare  various  finite  difference  schemes.  We  use  a  different  computational  domain  (Fig. 
4a)  from  the  one  used  in  [W].  Our  domain  is  physically  more  natural  and  computationally 
easily  manageable  for  triangulations.  This  is  in  contrast  with  finite  difference  methods 
which  would  meet  complications  in  using  our  domain.  The  initial  condition  is  described 
in  Fig.  4a.  It  corresponds  to  a  Mach  10  shock  making  60  degrees  with  the  bottom  wall. 
The  boundary  conditions  applied  are:  inflow  at  x  =  0;  outflow  at  x  =  3;  reflecting  at  the 
bottom;  and  the  exact  solution  of  a  Mach  10  shock  at  the  top.  In  Fig.  4c  we  show  the 
density  contour  at  T  =  0.2.  The  triangulation  used  is  drawn  in  Fig.  4b.  Again  ,  we  do  not 
try  to  align  the  triangles  with  the  shocks,  and  the  number  of  degrees  of  freedom  used  is 
close  to  the  middle  case  (Ax  =  Ay  =  1/60)  in  [W].  We  again  observe  a  result  comparable 
with  the  finite  difference  schemes,  except  for  the  above  mentioned  smearing  of  contact 
discontinuities. 


4.  Concluding  remarks.  We  have  discussed  the  generalization  of  the  RKDG  method 
to  two-dimensional  systems  of  conservation  laws,  using  the  Euler  equations  of  gas  dynamics 
as  an  example.  Numerical  tests  using  P1  elements  on  three  standard  examples  show 
comparable  results  with  nonoscillatory  finite  difference  schemes  and  indicate  good  potential 
of  our  finite  element  method  in  handling  geometry  and  boundary  conditions. 
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Fig.  3. a.  The  triangulation  ‘Ax  —  Ay  =  1/5’.  The  densities  below  hve  been 
computed  on  a  triangulation  eight  times  finer. 


Fig.  3.B.  Density  contours  computed  on  the  triangulation  ‘Ax  —  Ay  =  1/40’. 
No  treatment  of  the  singularity  at  the  comer  is  used. 


FlG.  3.C.  Density  contours  computed  on  the  triangulation  ‘Ax  =  Ay  —  1/40’. 
The  treatment  of  the  singularity  at  the  comer  suggested  in  (W ]  is  used. 


Fig.  4. a.  The  computational  domain  and  the  initial  conditions. 
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Fig.  4.B.  The  triangulation  ‘Ax  =  Ay  =  1/10’.  The  density  below  has  been 
computed  of  a  triangulation  six  times  finer. 


Fig.  4.c. 
1/60  ’. 


Density  contours  computed  with  the  triangulation  ‘Ax  =  Ay  = 
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